!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \brief routines that build the integrals of the Vxc potential calculated
!>      for the atomic code
! **************************************************************************************************
MODULE atom_xc

   USE atom_types,                      ONLY: atom_basis_type,&
                                              atom_type,&
                                              gth_pseudo,&
                                              opmat_type,&
                                              sgp_pseudo,&
                                              upf_pseudo
   USE atom_utils,                      ONLY: atom_core_density,&
                                              atom_density,&
                                              coulomb_potential_numeric,&
                                              integrate_grid,&
                                              numpot_matrix
   USE cp_files,                        ONLY: close_file,&
                                              get_unit_number,&
                                              open_file
   USE input_constants,                 ONLY: xc_none
   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
                                              section_vals_type,&
                                              section_vals_val_get
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: fourpi,&
                                              pi
   USE xc_atom,                         ONLY: xc_rho_set_atom_update
   USE xc_derivative_desc,              ONLY: &
        deriv_norm_drho, deriv_norm_drhoa, deriv_norm_drhob, deriv_rho, deriv_rhoa, deriv_rhob, &
        deriv_tau, deriv_tau_a, deriv_tau_b
   USE xc_derivative_set_types,         ONLY: xc_derivative_set_type,&
                                              xc_dset_create,&
                                              xc_dset_get_derivative,&
                                              xc_dset_release,&
                                              xc_dset_zero_all
   USE xc_derivative_types,             ONLY: xc_derivative_get,&
                                              xc_derivative_type
   USE xc_derivatives,                  ONLY: xc_functionals_eval,&
                                              xc_functionals_get_needs
   USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
   USE xc_rho_set_types,                ONLY: xc_rho_set_create,&
                                              xc_rho_set_release,&
                                              xc_rho_set_type
#include "./base/base_uses.f90"

   IMPLICIT NONE

   PRIVATE

   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_xc'

! ZMP added coulomb_potential_numeric
! ZMP added public subroutines calculate_atom_zmp, calculate_atom_ext_vxc
   PUBLIC :: calculate_atom_vxc_lda, calculate_atom_vxc_lsd, &
             calculate_atom_zmp, calculate_atom_ext_vxc
   PUBLIC :: atom_vxc_lda, atom_vxc_lsd, atom_dpot_lda

CONTAINS

! **************************************************************************************************
!> \brief ZMP subroutine for the calculation of the effective constraint
!>        potential in the atomic code.
!> \param ext_density  external electron density
!> \param atom         information about the atomic kind
!> \param lprint       save exchange-correlation potential to the file 'linear_potentials.dat'
!> \param xcmat        exchange-correlation potential matrix
!> \author D. Varsano [daniele.varsano@nano.cnr.it]
! **************************************************************************************************
   SUBROUTINE calculate_atom_zmp(ext_density, atom, lprint, xcmat)
      REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: ext_density
      TYPE(atom_type), INTENT(INOUT)                     :: atom
      LOGICAL                                            :: lprint
      TYPE(opmat_type), OPTIONAL, POINTER                :: xcmat

      CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_atom_zmp'

      INTEGER                                            :: extunit, handle, ir, nr, z
      REAL(KIND=dp)                                      :: int1, int2
      REAL(KIND=dp), DIMENSION(:), POINTER               :: deltarho, rho_dum, vxc, vxc1, vxc2
      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rho

      CALL timeset(routineN, handle)

      nr = atom%basis%grid%nr
      z = atom%z
      ALLOCATE (rho(nr, 1), vxc(nr), vxc1(nr), vxc2(nr), rho_dum(nr), deltarho(nr))
      CALL atom_density(rho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
      !Vxc1
      int1 = integrate_grid(ext_density, atom%basis%grid)
      int1 = fourpi*int1
      CALL coulomb_potential_numeric(vxc1, rho(:, 1), atom%basis%grid)

      vxc1 = -vxc1/z

      !Vxc2
      rho_dum = rho(:, 1)*int1/z
      deltarho = rho_dum - ext_density
      int2 = integrate_grid(deltarho, atom%basis%grid)
      CALL coulomb_potential_numeric(vxc2, deltarho, atom%basis%grid)

      vxc2 = vxc2*atom%lambda

      !Vxc
      vxc = vxc1 + vxc2

      atom%energy%exc = fourpi*integrate_grid(vxc, rho(:, 1), atom%basis%grid)
      atom%rho_diff_integral = fourpi*int2

      IF (lprint) THEN
         extunit = get_unit_number()
         CALL open_file(file_name="linear_potentials.dat", file_status="UNKNOWN", &
                        file_form="FORMATTED", file_action="WRITE", &
                        unit_number=extunit)

         WRITE (extunit, fmt='("#",T10,"R[bohr]",T36,"Rho[au]",T61,"DRho[au]",T86,"V_XC[au]",T111,"V_XC1[au]",T136,"V_XC2[au]")')
         DO ir = 1, nr
            WRITE (extunit, fmt='(T1,E24.15,T26,E24.15,T51,E24.15,T76,E24.15,T101,E24.15,T126,E24.15)') &
               atom%basis%grid%rad(ir), rho(ir, 1), deltarho(ir), vxc(ir), vxc1(ir), vxc2(ir)
         END DO
         CALL close_file(unit_number=extunit)
      END IF

      IF (PRESENT(xcmat)) CALL numpot_matrix(xcmat%op, vxc, atom%basis, 0)

      DEALLOCATE (rho, vxc, vxc1, vxc2, rho_dum, deltarho)

      CALL timestop(handle)

   END SUBROUTINE calculate_atom_zmp

! **************************************************************************************************
!> \brief ZMP subroutine for the scf external density from the external v_xc.
!> \param vxc     exchange-correlation potential
!> \param atom    information about the atomic kind
!> \param lprint  save exchange-correlation potential to the file 'linear_potentials.dat'
!> \param xcmat   exchange-correlation potential matrix
!> \author D. Varsano [daniele.varsano@nano.cnr.it]
! **************************************************************************************************
   SUBROUTINE calculate_atom_ext_vxc(vxc, atom, lprint, xcmat)
      REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: vxc
      TYPE(atom_type), INTENT(INOUT)                     :: atom
      LOGICAL, INTENT(in)                                :: lprint
      TYPE(opmat_type), INTENT(inout), OPTIONAL, POINTER :: xcmat

      CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_atom_ext_vxc'

      INTEGER                                            :: extunit, handle, ir, nr
      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rho

      CALL timeset(routineN, handle)
      nr = atom%basis%grid%nr

      ALLOCATE (rho(nr, 1))

      CALL atom_density(rho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")

      IF (lprint) THEN
         extunit = get_unit_number()
         CALL open_file(file_name="linear_potentials.dat", file_status="UNKNOWN", &
                        file_form="FORMATTED", file_action="WRITE", &
                        unit_number=extunit)

         WRITE (extunit, fmt='("#",T10,"R[bohr]",T36,"Rho[au]",T61,"V_XC[au]")')
         DO ir = 1, nr
            WRITE (extunit, fmt='(T1,E24.15,T26,E24.15,T51,E24.15)') &
               atom%basis%grid%rad(ir), rho(ir, 1), vxc(ir)
         END DO
         CALL close_file(unit_number=extunit)
      END IF

      atom%energy%exc = fourpi*integrate_grid(vxc, rho(:, 1), atom%basis%grid)
      IF (PRESENT(xcmat)) CALL numpot_matrix(xcmat%op, vxc, atom%basis, 0)

      DEALLOCATE (rho)

      CALL timestop(handle)

   END SUBROUTINE calculate_atom_ext_vxc

! **************************************************************************************************
!> \brief Calculate atomic exchange-correlation potential in spin-restricted case.
!> \param xcmat       exchange-correlation potential expressed in the atomic basis set
!> \param atom        information about the atomic kind
!> \param xc_section  XC input section
!> \par History
!>    * 03.2010 extension of GTH pseudo-potential definition [Juerg Hutter]
!>    * 02.2010 renamed to calculate_atom_vxc_lda() [Juerg Hutter]
!>    * 08.2008 created [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE calculate_atom_vxc_lda(xcmat, atom, xc_section)
      TYPE(opmat_type), POINTER                          :: xcmat
      TYPE(atom_type), INTENT(INOUT)                     :: atom
      TYPE(section_vals_type), POINTER                   :: xc_section

      CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_atom_vxc_lda'

      INTEGER                                            :: deriv_order, handle, i, l, myfun, n1, &
                                                            n2, n3, nr, nspins, unit_nr
      INTEGER, DIMENSION(2, 3)                           :: bounds
      LOGICAL                                            :: lsd, nlcc
      REAL(KIND=dp)                                      :: density_cut, gradient_cut, tau_cut
      REAL(KIND=dp), DIMENSION(:), POINTER               :: exc, vxc
      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: drho, lap, rho, tau
      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: taumat, xcpot
      TYPE(section_vals_type), POINTER                   :: xc_fun_section
      TYPE(xc_derivative_set_type)                       :: deriv_set
      TYPE(xc_derivative_type), POINTER                  :: deriv
      TYPE(xc_rho_cflags_type)                           :: needs
      TYPE(xc_rho_set_type)                              :: rho_set

      CALL timeset(routineN, handle)

      nlcc = .FALSE.
      IF (atom%potential%ppot_type == gth_pseudo) THEN
         nlcc = atom%potential%gth_pot%nlcc
      ELSE IF (atom%potential%ppot_type == upf_pseudo) THEN
         nlcc = atom%potential%upf_pot%core_correction
      ELSE IF (atom%potential%ppot_type == sgp_pseudo) THEN
         nlcc = atom%potential%sgp_pot%has_nlcc
      END IF

      IF (ASSOCIATED(xc_section)) THEN
         xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
         CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)

         IF (myfun == xc_none) THEN
            atom%energy%exc = 0._dp
         ELSE
            CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
            CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
            CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)

            lsd = .FALSE.
            nspins = 1
            needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.FALSE.)

            ! Prepare the structures needed to calculate and store the xc derivatives

            ! Array dimension: here anly one dimensional arrays are used,
            ! i.e. only the first column of deriv_data is read.
            ! The other to dimensions  are set to size equal 1
            nr = atom%basis%grid%nr
            bounds(1:2, 1:3) = 1
            bounds(2, 1) = nr

            ! create a place where to put the derivatives
            CALL xc_dset_create(deriv_set, local_bounds=bounds)
            ! create the place where to store the argument for the functionals
            CALL xc_rho_set_create(rho_set, bounds, rho_cutoff=density_cut, &
                                   drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
            ! allocate the required 3d arrays where to store rho and drho
            CALL xc_rho_set_atom_update(rho_set, needs, nspins, bounds)

            NULLIFY (rho, drho, tau)
            IF (needs%rho) THEN
               ALLOCATE (rho(nr, 1))
               CALL atom_density(rho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
               IF (nlcc) THEN
                  CALL atom_core_density(rho(:, 1), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
               END IF
            END IF
            IF (needs%norm_drho) THEN
               ALLOCATE (drho(nr, 1))
               CALL atom_density(drho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="DER")
               IF (nlcc) THEN
                  CALL atom_core_density(drho(:, 1), atom%potential, typ="DER", rr=atom%basis%grid%rad)
               END IF
            END IF
            IF (needs%tau) THEN
               ALLOCATE (tau(nr, 1))
               CALL atom_density(tau(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, &
                                 typ="KIN", rr=atom%basis%grid%rad2)
            END IF
            IF (needs%laplace_rho) THEN
               ALLOCATE (lap(nr, 1))
               CALL atom_density(lap(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, &
                                 typ="LAP", rr=atom%basis%grid%rad2)
               IF (nlcc) THEN
                  CALL atom_core_density(lap(:, 1), atom%potential, typ="LAP", rr=atom%basis%grid%rad)
               END IF
            END IF

            CALL fill_rho_set(rho_set, nspins, needs, rho, drho, tau, lap, nr)

            CALL xc_dset_zero_all(deriv_set)

            deriv_order = 1
            CALL xc_functionals_eval(xc_fun_section, lsd=lsd, rho_set=rho_set, deriv_set=deriv_set, &
                                     deriv_order=deriv_order)

            ! Integration to get the matrix elements and energy
            deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], allocate_deriv=.FALSE.)
            CALL xc_derivative_get(deriv, deriv_data=xcpot)
            atom%energy%exc = fourpi*integrate_grid(xcpot(:, 1, 1), atom%basis%grid)

            ! dump grid density and xcpot (xc energy?)
            IF (.FALSE.) THEN
               CALL open_file(unit_number=unit_nr, file_status="UNKNOWN", file_action="WRITE", file_name="atom.dat")
               DO i = 1, SIZE(xcpot(:, 1, 1))
                  WRITE (unit_nr, *) atom%basis%grid%rad(i), rho(i, 1), xcpot(i, 1, 1)
               END DO
               CALL close_file(unit_nr)
            END IF

            IF (needs%rho) THEN
               deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], allocate_deriv=.FALSE.)
               CALL xc_derivative_get(deriv, deriv_data=xcpot)
               CALL numpot_matrix(xcmat%op, xcpot(:, 1, 1), atom%basis, 0)
               IF (.FALSE.) THEN
                  CALL open_file(unit_number=unit_nr, file_status="UNKNOWN", file_action="WRITE", file_name="atompot.dat")
                  DO i = 1, SIZE(xcpot(:, 1, 1))
                     WRITE (unit_nr, *) atom%basis%grid%rad(i), rho(i, 1), xcpot(i, 1, 1)
                  END DO
                  CALL close_file(unit_nr)
               END IF
               DEALLOCATE (rho)
            END IF
            IF (needs%norm_drho) THEN
               deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], allocate_deriv=.FALSE.)
               CALL xc_derivative_get(deriv, deriv_data=xcpot)
               CALL numpot_matrix(xcmat%op, xcpot(:, 1, 1), atom%basis, 1)
               IF (.FALSE.) THEN
                  CALL open_file(unit_number=unit_nr, file_status="UNKNOWN", file_action="WRITE", file_name="atomdpot.dat")
                  DO i = 1, SIZE(xcpot(:, 1, 1))
                     WRITE (unit_nr, *) atom%basis%grid%rad(i), drho(i, 1), xcpot(i, 1, 1)
                  END DO
                  CALL close_file(unit_nr)
               END IF
               DEALLOCATE (drho)
            END IF
            IF (needs%tau) THEN
               deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], allocate_deriv=.FALSE.)
               CALL xc_derivative_get(deriv, deriv_data=xcpot)
               n1 = SIZE(xcmat%op, 1)
               n2 = SIZE(xcmat%op, 2)
               n3 = SIZE(xcmat%op, 3)
               ALLOCATE (taumat(n1, n2, 0:n3 - 1))
               taumat = 0._dp

               xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
               CALL numpot_matrix(xcmat%op, xcpot(:, 1, 1), atom%basis, 2)
               xcpot(:, 1, 1) = xcpot(:, 1, 1)/atom%basis%grid%rad2(:)
               CALL numpot_matrix(taumat, xcpot(:, 1, 1), atom%basis, 0)
               DO l = 0, 3
                  xcmat%op(:, :, l) = xcmat%op(:, :, l) + REAL(l*(l + 1), dp)*taumat(:, :, l)
               END DO

               DEALLOCATE (tau)
               DEALLOCATE (taumat)
            END IF

            ! assume lap is not really needed
            IF (needs%laplace_rho) THEN
               DEALLOCATE (lap)
            END IF

            ! Release the xc structure used to store the xc derivatives
            CALL xc_dset_release(deriv_set)
            CALL xc_rho_set_release(rho_set)

         END IF !xc_none

      ELSE

         ! we don't have an xc_section, use a default setup
         nr = atom%basis%grid%nr
         ALLOCATE (rho(nr, 1), exc(nr), vxc(nr))

         CALL atom_density(rho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
         IF (nlcc) THEN
            CALL atom_core_density(rho(:, 1), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
         END IF
         CALL lda_pade(rho(:, 1), exc, vxc)

         atom%energy%exc = fourpi*integrate_grid(exc, atom%basis%grid)
         CALL numpot_matrix(xcmat%op, vxc, atom%basis, 0)

         DEALLOCATE (rho, exc, vxc)

      END IF

      CALL timestop(handle)

   END SUBROUTINE calculate_atom_vxc_lda

! **************************************************************************************************
!> \brief Calculate atomic exchange-correlation potential in spin-polarized case.
!> \param xcmata      spin-alpha exchange-correlation potential matrix
!> \param xcmatb      spin-beta exchange-correlation potential matrix
!> \param atom        information about the atomic kind
!> \param xc_section  XC input section
!> \par History
!>    * 02.2010 created [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE calculate_atom_vxc_lsd(xcmata, xcmatb, atom, xc_section)
      TYPE(opmat_type), POINTER                          :: xcmata, xcmatb
      TYPE(atom_type), INTENT(INOUT)                     :: atom
      TYPE(section_vals_type), POINTER                   :: xc_section

      CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_atom_vxc_lsd'

      INTEGER                                            :: deriv_order, handle, l, myfun, n1, n2, &
                                                            n3, nr, nspins
      INTEGER, DIMENSION(2, 3)                           :: bounds
      LOGICAL                                            :: lsd, nlcc
      REAL(KIND=dp)                                      :: density_cut, gradient_cut, tau_cut
      REAL(KIND=dp), DIMENSION(:), POINTER               :: exc, vxca, vxcb, xfun
      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: drho, lap, rho, tau
      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: taumat, xcpot
      TYPE(section_vals_type), POINTER                   :: xc_fun_section
      TYPE(xc_derivative_set_type)                       :: deriv_set
      TYPE(xc_derivative_type), POINTER                  :: deriv
      TYPE(xc_rho_cflags_type)                           :: needs
      TYPE(xc_rho_set_type)                              :: rho_set

      CALL timeset(routineN, handle)

      nlcc = .FALSE.
      IF (atom%potential%ppot_type == gth_pseudo) THEN
         nlcc = atom%potential%gth_pot%nlcc
      ELSE IF (atom%potential%ppot_type == upf_pseudo) THEN
         nlcc = atom%potential%upf_pot%core_correction
         CPASSERT(.NOT. nlcc)
      ELSE IF (atom%potential%ppot_type == sgp_pseudo) THEN
         nlcc = atom%potential%sgp_pot%has_nlcc
         CPASSERT(.NOT. nlcc)
      END IF

      IF (ASSOCIATED(xc_section)) THEN
         IF (nlcc) THEN
            nr = atom%basis%grid%nr
            ALLOCATE (xfun(nr))
         END IF

         xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
         CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)

         IF (myfun == xc_none) THEN
            atom%energy%exc = 0._dp
         ELSE
            CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
            CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
            CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)

            lsd = .TRUE.
            nspins = 2
            needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.FALSE.)

            ! Prepare the structures needed to calculate and store the xc derivatives

            ! Array dimension: here anly one dimensional arrays are used,
            ! i.e. only the first column of deriv_data is read.
            ! The other to dimensions  are set to size equal 1
            nr = atom%basis%grid%nr
            bounds(1:2, 1:3) = 1
            bounds(2, 1) = nr

            ! create a place where to put the derivatives
            CALL xc_dset_create(deriv_set, local_bounds=bounds)
            ! create the place where to store the argument for the functionals
            CALL xc_rho_set_create(rho_set, bounds, rho_cutoff=density_cut, &
                                   drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
            ! allocate the required 3d arrays where to store rho and drho
            CALL xc_rho_set_atom_update(rho_set, needs, nspins, bounds)

            NULLIFY (rho, drho, tau)
            IF (needs%rho_spin) THEN
               ALLOCATE (rho(nr, 2))
               CALL atom_density(rho(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, typ="RHO")
               CALL atom_density(rho(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, typ="RHO")
               IF (nlcc) THEN
                  xfun = 0.0_dp
                  CALL atom_core_density(xfun(:), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
                  rho(:, 1) = rho(:, 1) + 0.5_dp*xfun(:)
                  rho(:, 2) = rho(:, 2) + 0.5_dp*xfun(:)
               END IF
            END IF
            IF (needs%norm_drho_spin) THEN
               ALLOCATE (drho(nr, 2))
               CALL atom_density(drho(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, typ="DER")
               CALL atom_density(drho(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, typ="DER")
               IF (nlcc) THEN
                  xfun = 0.0_dp
                  CALL atom_core_density(xfun(:), atom%potential, typ="DER", rr=atom%basis%grid%rad)
                  drho(:, 1) = drho(:, 1) + 0.5_dp*xfun(:)
                  drho(:, 2) = drho(:, 2) + 0.5_dp*xfun(:)
               END IF
            END IF
            IF (needs%tau_spin) THEN
               ALLOCATE (tau(nr, 2))
               CALL atom_density(tau(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, &
                                 typ="KIN", rr=atom%basis%grid%rad2)
               CALL atom_density(tau(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, &
                                 typ="KIN", rr=atom%basis%grid%rad2)
            END IF
            IF (needs%laplace_rho_spin) THEN
               ALLOCATE (lap(nr, 2))
               CALL atom_density(lap(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, &
                                 typ="LAP", rr=atom%basis%grid%rad2)
               CALL atom_density(lap(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, &
                                 typ="LAP", rr=atom%basis%grid%rad2)
               IF (nlcc) THEN
                  xfun = 0.0_dp
                  CALL atom_core_density(xfun(:), atom%potential, typ="LAP", rr=atom%basis%grid%rad)
                  lap(:, 1) = lap(:, 1) + 0.5_dp*xfun(:)
                  lap(:, 2) = lap(:, 2) + 0.5_dp*xfun(:)
               END IF
            END IF

            CALL fill_rho_set(rho_set, nspins, needs, rho, drho, tau, lap, nr)

            CALL xc_dset_zero_all(deriv_set)

            deriv_order = 1
            CALL xc_functionals_eval(xc_fun_section, lsd=lsd, rho_set=rho_set, deriv_set=deriv_set, &
                                     deriv_order=deriv_order)

            ! Integration to get the matrix elements and energy
            deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], allocate_deriv=.FALSE.)
            CALL xc_derivative_get(deriv, deriv_data=xcpot)
            atom%energy%exc = fourpi*integrate_grid(xcpot(:, 1, 1), atom%basis%grid)

            IF (needs%rho_spin) THEN
               deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], allocate_deriv=.FALSE.)
               CALL xc_derivative_get(deriv, deriv_data=xcpot)
               CALL numpot_matrix(xcmata%op, xcpot(:, 1, 1), atom%basis, 0)
               deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], allocate_deriv=.FALSE.)
               CALL xc_derivative_get(deriv, deriv_data=xcpot)
               CALL numpot_matrix(xcmatb%op, xcpot(:, 1, 1), atom%basis, 0)
               DEALLOCATE (rho)
            END IF
            IF (needs%norm_drho_spin) THEN
               ! drhoa
               NULLIFY (deriv)
               deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], allocate_deriv=.FALSE.)
               CALL xc_derivative_get(deriv, deriv_data=xcpot)
               CALL numpot_matrix(xcmata%op, xcpot(:, 1, 1), atom%basis, 1)
               ! drhob
               NULLIFY (deriv)
               deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], allocate_deriv=.FALSE.)
               CALL xc_derivative_get(deriv, deriv_data=xcpot)
               CALL numpot_matrix(xcmatb%op, xcpot(:, 1, 1), atom%basis, 1)
               ! Cross Terms
               NULLIFY (deriv)
               deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
               IF (ASSOCIATED(deriv)) THEN
                  CALL xc_derivative_get(deriv, deriv_data=xcpot)
                  CALL numpot_matrix(xcmata%op, xcpot(:, 1, 1), atom%basis, 1)
                  CALL numpot_matrix(xcmatb%op, xcpot(:, 1, 1), atom%basis, 1)
               END IF
               DEALLOCATE (drho)
            END IF
            IF (needs%tau_spin) THEN
               n1 = SIZE(xcmata%op, 1)
               n2 = SIZE(xcmata%op, 2)
               n3 = SIZE(xcmata%op, 3)
               ALLOCATE (taumat(n1, n2, 0:n3 - 1))

               deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a], allocate_deriv=.FALSE.)
               CALL xc_derivative_get(deriv, deriv_data=xcpot)
               taumat = 0._dp
               xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
               CALL numpot_matrix(xcmata%op, xcpot(:, 1, 1), atom%basis, 2)
               xcpot(:, 1, 1) = xcpot(:, 1, 1)/atom%basis%grid%rad2(:)
               CALL numpot_matrix(taumat, xcpot(:, 1, 1), atom%basis, 0)
               DO l = 0, 3
                  xcmata%op(:, :, l) = xcmata%op(:, :, l) + REAL(l*(l + 1), dp)*taumat(:, :, l)
               END DO

               deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b], allocate_deriv=.FALSE.)
               CALL xc_derivative_get(deriv, deriv_data=xcpot)
               taumat = 0._dp
               xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
               CALL numpot_matrix(xcmatb%op, xcpot(:, 1, 1), atom%basis, 2)
               xcpot(:, 1, 1) = xcpot(:, 1, 1)/atom%basis%grid%rad2(:)
               CALL numpot_matrix(taumat, xcpot(:, 1, 1), atom%basis, 0)
               DO l = 0, 3
                  xcmatb%op(:, :, l) = xcmatb%op(:, :, l) + REAL(l*(l + 1), dp)*taumat(:, :, l)
               END DO

               DEALLOCATE (tau)
               DEALLOCATE (taumat)
            END IF

            ! assume lap is not really needed
            IF (needs%laplace_rho_spin) THEN
               DEALLOCATE (lap)
            END IF

            ! Release the xc structure used to store the xc derivatives
            CALL xc_dset_release(deriv_set)
            CALL xc_rho_set_release(rho_set)

         END IF !xc_none

         IF (nlcc) DEALLOCATE (xfun)

      ELSE

         ! we don't have an xc_section, use a default setup
         nr = atom%basis%grid%nr
         ALLOCATE (rho(nr, 2), exc(nr), vxca(nr), vxcb(nr))
         IF (nlcc) ALLOCATE (xfun(nr))

         CALL atom_density(rho(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, typ="RHO")
         CALL atom_density(rho(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, typ="RHO")
         IF (nlcc) THEN
            xfun(:) = 0.0_dp
            CALL atom_core_density(xfun(:), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
            rho(:, 1) = rho(:, 1) + 0.5_dp*xfun(:)
            rho(:, 2) = rho(:, 2) + 0.5_dp*xfun(:)
         END IF
         CALL lsd_pade(rho(:, 1), rho(:, 2), exc, vxca, vxcb)

         atom%energy%exc = fourpi*integrate_grid(exc, atom%basis%grid)
         CALL numpot_matrix(xcmata%op, vxca, atom%basis, 0)
         CALL numpot_matrix(xcmatb%op, vxcb, atom%basis, 0)

         DEALLOCATE (rho, exc, vxca, vxcb)
         IF (nlcc) DEALLOCATE (xfun)

      END IF

      CALL timestop(handle)

   END SUBROUTINE calculate_atom_vxc_lsd

! **************************************************************************************************
!> \brief Alternative subroutine that calculates atomic exchange-correlation potential
!>        in spin-restricted case.
!> \param basis       atomic basis set
!> \param pmat        density matrix
!> \param maxl        maximum angular momentum
!> \param xc_section  XC input section
!> \param fexc        exchange-correlation energy
!> \param xcmat       exchange-correlation potential matrix
!> \par History
!>    * 07.2016 ADMM analysis [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE atom_vxc_lda(basis, pmat, maxl, xc_section, fexc, xcmat)
      TYPE(atom_basis_type), POINTER                     :: basis
      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: pmat
      INTEGER, INTENT(IN)                                :: maxl
      TYPE(section_vals_type), POINTER                   :: xc_section
      REAL(KIND=dp), INTENT(OUT)                         :: fexc
      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: xcmat

      CHARACTER(LEN=*), PARAMETER                        :: routineN = 'atom_vxc_lda'

      INTEGER                                            :: deriv_order, handle, l, myfun, n1, n2, &
                                                            n3, nr, nspins
      INTEGER, DIMENSION(2, 3)                           :: bounds
      LOGICAL                                            :: lsd
      REAL(KIND=dp)                                      :: density_cut, gradient_cut, tau_cut
      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: drho, lap, rho, tau
      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: taumat, xcpot
      TYPE(section_vals_type), POINTER                   :: xc_fun_section
      TYPE(xc_derivative_set_type)                       :: deriv_set
      TYPE(xc_derivative_type), POINTER                  :: deriv
      TYPE(xc_rho_cflags_type)                           :: needs
      TYPE(xc_rho_set_type)                              :: rho_set

      CALL timeset(routineN, handle)

      CPASSERT(ASSOCIATED(xc_section))

      xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
      CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)

      IF (myfun == xc_none) THEN
         fexc = 0._dp
      ELSE
         CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
         CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
         CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)

         lsd = .FALSE.
         nspins = 1
         needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.FALSE.)

         ! Prepare the structures needed to calculate and store the xc derivatives

         ! Array dimension: here anly one dimensional arrays are used,
         ! i.e. only the first column of deriv_data is read.
         ! The other to dimensions  are set to size equal 1
         nr = basis%grid%nr
         bounds(1:2, 1:3) = 1
         bounds(2, 1) = nr

         ! create a place where to put the derivatives
         CALL xc_dset_create(deriv_set, local_bounds=bounds)
         ! create the place where to store the argument for the functionals
         CALL xc_rho_set_create(rho_set, bounds, rho_cutoff=density_cut, &
                                drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
         ! allocate the required 3d arrays where to store rho and drho
         CALL xc_rho_set_atom_update(rho_set, needs, nspins, bounds)

         NULLIFY (rho, drho, tau)
         IF (needs%rho) THEN
            ALLOCATE (rho(nr, 1))
            CALL atom_density(rho(:, 1), pmat, basis, maxl, typ="RHO")
         END IF
         IF (needs%norm_drho) THEN
            ALLOCATE (drho(nr, 1))
            CALL atom_density(drho(:, 1), pmat, basis, maxl, typ="DER")
         END IF
         IF (needs%tau) THEN
            ALLOCATE (tau(nr, 1))
            CALL atom_density(tau(:, 1), pmat, basis, maxl, typ="KIN", rr=basis%grid%rad2)
         END IF
         IF (needs%laplace_rho) THEN
            ALLOCATE (lap(nr, 1))
            CALL atom_density(lap(:, 1), pmat, basis, maxl, typ="LAP", rr=basis%grid%rad2)
         END IF

         CALL fill_rho_set(rho_set, nspins, needs, rho, drho, tau, lap, nr)

         CALL xc_dset_zero_all(deriv_set)

         deriv_order = 1
         CALL xc_functionals_eval(xc_fun_section, lsd=lsd, rho_set=rho_set, deriv_set=deriv_set, &
                                  deriv_order=deriv_order)

         ! Integration to get the matrix elements and energy
         deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], allocate_deriv=.FALSE.)
         CALL xc_derivative_get(deriv, deriv_data=xcpot)
         fexc = fourpi*integrate_grid(xcpot(:, 1, 1), basis%grid)

         IF (needs%rho) THEN
            deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], allocate_deriv=.FALSE.)
            CALL xc_derivative_get(deriv, deriv_data=xcpot)
            CALL numpot_matrix(xcmat, xcpot(:, 1, 1), basis, 0)
            DEALLOCATE (rho)
         END IF
         IF (needs%norm_drho) THEN
            deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], allocate_deriv=.FALSE.)
            CALL xc_derivative_get(deriv, deriv_data=xcpot)
            CALL numpot_matrix(xcmat, xcpot(:, 1, 1), basis, 1)
            DEALLOCATE (drho)
         END IF
         IF (needs%tau) THEN
            deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], allocate_deriv=.FALSE.)
            CALL xc_derivative_get(deriv, deriv_data=xcpot)
            n1 = SIZE(xcmat, 1)
            n2 = SIZE(xcmat, 2)
            n3 = SIZE(xcmat, 3)
            ALLOCATE (taumat(n1, n2, 0:n3 - 1))
            taumat = 0._dp

            xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
            CALL numpot_matrix(xcmat, xcpot(:, 1, 1), basis, 2)
            xcpot(:, 1, 1) = xcpot(:, 1, 1)/basis%grid%rad2(:)
            CALL numpot_matrix(taumat, xcpot(:, 1, 1), basis, 0)
            DO l = 0, 3
               xcmat(:, :, l) = xcmat(:, :, l) + REAL(l*(l + 1), dp)*taumat(:, :, l)
            END DO

            DEALLOCATE (tau)
            DEALLOCATE (taumat)
         END IF

         ! we assume lap is not really needed
         IF (needs%laplace_rho) THEN
            DEALLOCATE (lap)
         END IF

         ! Release the xc structure used to store the xc derivatives
         CALL xc_dset_release(deriv_set)
         CALL xc_rho_set_release(rho_set)

      END IF !xc_none

      CALL timestop(handle)

   END SUBROUTINE atom_vxc_lda

! **************************************************************************************************
!> \brief Alternative subroutine that calculates atomic exchange-correlation potential
!>        in spin-polarized case.
!> \param basis       atomic basis set
!> \param pmata       spin-alpha density matrix
!> \param pmatb       spin-beta density matrix
!> \param maxl        maximum angular momentum
!> \param xc_section  XC input section
!> \param fexc        exchange-correlation energy
!> \param xcmata      spin-alpha exchange-correlation potential matrix
!> \param xcmatb      spin-beta exchange-correlation potential matrix
!> \par History
!>    * 07.2016 ADMM analysis [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE atom_vxc_lsd(basis, pmata, pmatb, maxl, xc_section, fexc, xcmata, xcmatb)
      TYPE(atom_basis_type), POINTER                     :: basis
      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: pmata, pmatb
      INTEGER, INTENT(IN)                                :: maxl
      TYPE(section_vals_type), POINTER                   :: xc_section
      REAL(KIND=dp), INTENT(OUT)                         :: fexc
      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: xcmata, xcmatb

      CHARACTER(LEN=*), PARAMETER                        :: routineN = 'atom_vxc_lsd'

      INTEGER                                            :: deriv_order, handle, l, myfun, n1, n2, &
                                                            n3, nr, nspins
      INTEGER, DIMENSION(2, 3)                           :: bounds
      LOGICAL                                            :: lsd
      REAL(KIND=dp)                                      :: density_cut, gradient_cut, tau_cut
      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: drho, lap, rho, tau
      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: taumat, xcpot
      TYPE(section_vals_type), POINTER                   :: xc_fun_section
      TYPE(xc_derivative_set_type)                       :: deriv_set
      TYPE(xc_derivative_type), POINTER                  :: deriv
      TYPE(xc_rho_cflags_type)                           :: needs
      TYPE(xc_rho_set_type)                              :: rho_set

      CALL timeset(routineN, handle)

      CPASSERT(ASSOCIATED(xc_section))

      xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
      CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)

      IF (myfun == xc_none) THEN
         fexc = 0._dp
      ELSE
         CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
         CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
         CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)

         lsd = .TRUE.
         nspins = 2
         needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.FALSE.)

         ! Prepare the structures needed to calculate and store the xc derivatives

         ! Array dimension: here anly one dimensional arrays are used,
         ! i.e. only the first column of deriv_data is read.
         ! The other to dimensions  are set to size equal 1
         nr = basis%grid%nr
         bounds(1:2, 1:3) = 1
         bounds(2, 1) = nr

         ! create a place where to put the derivatives
         CALL xc_dset_create(deriv_set, local_bounds=bounds)
         ! create the place where to store the argument for the functionals
         CALL xc_rho_set_create(rho_set, bounds, rho_cutoff=density_cut, &
                                drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
         ! allocate the required 3d arrays where to store rho and drho
         CALL xc_rho_set_atom_update(rho_set, needs, nspins, bounds)

         NULLIFY (rho, drho, tau)
         IF (needs%rho_spin) THEN
            ALLOCATE (rho(nr, 2))
            CALL atom_density(rho(:, 1), pmata, basis, maxl, typ="RHO")
            CALL atom_density(rho(:, 2), pmatb, basis, maxl, typ="RHO")
         END IF
         IF (needs%norm_drho_spin) THEN
            ALLOCATE (drho(nr, 2))
            CALL atom_density(drho(:, 1), pmata, basis, maxl, typ="DER")
            CALL atom_density(drho(:, 2), pmatb, basis, maxl, typ="DER")
         END IF
         IF (needs%tau_spin) THEN
            ALLOCATE (tau(nr, 2))
            CALL atom_density(tau(:, 1), pmata, basis, maxl, typ="KIN", rr=basis%grid%rad2)
            CALL atom_density(tau(:, 2), pmatb, basis, maxl, typ="KIN", rr=basis%grid%rad2)
         END IF
         IF (needs%laplace_rho_spin) THEN
            ALLOCATE (lap(nr, 2))
            CALL atom_density(lap(:, 1), pmata, basis, maxl, typ="LAP", rr=basis%grid%rad2)
            CALL atom_density(lap(:, 2), pmatb, basis, maxl, typ="LAP", rr=basis%grid%rad2)
         END IF

         CALL fill_rho_set(rho_set, nspins, needs, rho, drho, tau, lap, nr)

         CALL xc_dset_zero_all(deriv_set)

         deriv_order = 1
         CALL xc_functionals_eval(xc_fun_section, lsd=lsd, rho_set=rho_set, deriv_set=deriv_set, &
                                  deriv_order=deriv_order)

         ! Integration to get the matrix elements and energy
         deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], allocate_deriv=.FALSE.)
         CALL xc_derivative_get(deriv, deriv_data=xcpot)
         fexc = fourpi*integrate_grid(xcpot(:, 1, 1), basis%grid)

         IF (needs%rho_spin) THEN
            deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], allocate_deriv=.FALSE.)
            CALL xc_derivative_get(deriv, deriv_data=xcpot)
            CALL numpot_matrix(xcmata, xcpot(:, 1, 1), basis, 0)
            deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], allocate_deriv=.FALSE.)
            CALL xc_derivative_get(deriv, deriv_data=xcpot)
            CALL numpot_matrix(xcmatb, xcpot(:, 1, 1), basis, 0)
            DEALLOCATE (rho)
         END IF
         IF (needs%norm_drho_spin) THEN
            ! drhoa
            NULLIFY (deriv)
            deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], allocate_deriv=.FALSE.)
            CALL xc_derivative_get(deriv, deriv_data=xcpot)
            CALL numpot_matrix(xcmata, xcpot(:, 1, 1), basis, 1)
            ! drhob
            NULLIFY (deriv)
            deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], allocate_deriv=.FALSE.)
            CALL xc_derivative_get(deriv, deriv_data=xcpot)
            CALL numpot_matrix(xcmatb, xcpot(:, 1, 1), basis, 1)
            ! Cross Terms
            NULLIFY (deriv)
            deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
            IF (ASSOCIATED(deriv)) THEN
               CALL xc_derivative_get(deriv, deriv_data=xcpot)
               CALL numpot_matrix(xcmata, xcpot(:, 1, 1), basis, 1)
               CALL numpot_matrix(xcmatb, xcpot(:, 1, 1), basis, 1)
            END IF
            DEALLOCATE (drho)
         END IF
         IF (needs%tau_spin) THEN
            n1 = SIZE(xcmata, 1)
            n2 = SIZE(xcmata, 2)
            n3 = SIZE(xcmata, 3)
            ALLOCATE (taumat(n1, n2, 0:n3 - 1))

            deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a], allocate_deriv=.FALSE.)
            CALL xc_derivative_get(deriv, deriv_data=xcpot)
            taumat = 0._dp
            xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
            CALL numpot_matrix(xcmata, xcpot(:, 1, 1), basis, 2)
            xcpot(:, 1, 1) = xcpot(:, 1, 1)/basis%grid%rad2(:)
            CALL numpot_matrix(taumat, xcpot(:, 1, 1), basis, 0)
            DO l = 0, 3
               xcmata(:, :, l) = xcmata(:, :, l) + REAL(l*(l + 1), dp)*taumat(:, :, l)
            END DO

            deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b], allocate_deriv=.FALSE.)
            CALL xc_derivative_get(deriv, deriv_data=xcpot)
            taumat = 0._dp
            xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
            CALL numpot_matrix(xcmatb, xcpot(:, 1, 1), basis, 2)
            xcpot(:, 1, 1) = xcpot(:, 1, 1)/basis%grid%rad2(:)
            CALL numpot_matrix(taumat, xcpot(:, 1, 1), basis, 0)
            DO l = 0, 3
               xcmatb(:, :, l) = xcmatb(:, :, l) + REAL(l*(l + 1), dp)*taumat(:, :, l)
            END DO

            DEALLOCATE (tau)
            DEALLOCATE (taumat)
         END IF

         ! Release the xc structure used to store the xc derivatives
         CALL xc_dset_release(deriv_set)
         CALL xc_rho_set_release(rho_set)

      END IF !xc_none

      CALL timestop(handle)

   END SUBROUTINE atom_vxc_lsd

! **************************************************************************************************
!> \brief Estimate the ADMM exchange energy correction term using a model exchange functional.
!> \param basis0      reference atomic basis set
!> \param pmat0       reference density matrix
!> \param basis1      ADMM basis set
!> \param pmat1       ADMM density matrix
!> \param maxl        maxumum angular momentum
!> \param functional  name of the model exchange functional:
!>                    "LINX" -- linear extrapolation of the Slater exchange potential [?]
!> \param dfexc       exchange-correlation energy difference
!> \param linxpar     LINX parameters
!> \par History
!>    * 07.2016 ADMM analysis [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE atom_dpot_lda(basis0, pmat0, basis1, pmat1, maxl, functional, dfexc, linxpar)
      TYPE(atom_basis_type), POINTER                     :: basis0
      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: pmat0
      TYPE(atom_basis_type), POINTER                     :: basis1
      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: pmat1
      INTEGER, INTENT(IN)                                :: maxl
      CHARACTER(LEN=*), INTENT(IN)                       :: functional
      REAL(KIND=dp), INTENT(OUT)                         :: dfexc
      REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL  :: linxpar

      CHARACTER(LEN=*), PARAMETER                        :: routineN = 'atom_dpot_lda'
      REAL(KIND=dp), PARAMETER :: alx = 0.20_dp, blx = -0.06_dp, &
         fs = -0.75_dp*(3._dp/pi)**(1._dp/3._dp)

      INTEGER                                            :: handle, ir, nr
      REAL(KIND=dp)                                      :: a, b, fx
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: delta, drho0, drho1, pot0, pot1, rho0, &
                                                            rho1

      CALL timeset(routineN, handle)

      nr = basis0%grid%nr

      ALLOCATE (rho0(nr), drho0(nr))
      CALL atom_density(rho0, pmat0, basis0, maxl, typ="RHO")
      CALL atom_density(drho0, pmat0, basis0, maxl, typ="DER")
      !
      ALLOCATE (rho1(nr), drho1(nr))
      CALL atom_density(rho1, pmat1, basis1, maxl, typ="RHO")
      ! CONSIDER TO REMOVE [?]: drho1 is calculated but it is not used.
      CALL atom_density(drho1, pmat1, basis1, maxl, typ="DER")
      !
      ALLOCATE (delta(nr))
      fx = 4.0_dp/3.0_dp
      delta(1:nr) = fs*(rho0(1:nr)**fx - rho1(1:nr)**fx)

      SELECT CASE (TRIM(functional))
      CASE ("LINX")
         IF (PRESENT(linxpar)) THEN
            a = linxpar(1)
            b = linxpar(2)
         ELSE
            a = alx
            b = blx
         END IF
         ALLOCATE (pot0(nr), pot1(nr))
         DO ir = 1, nr
            IF (rho0(ir) > 1.e-12_dp) THEN
               pot0(ir) = 0.5_dp*drho0(ir)/(3._dp*pi*pi*rho0(ir)**fx)
            ELSE
               pot0(ir) = 0._dp
            END IF
         END DO
         pot1(1:nr) = 1._dp + (a*pot0(1:nr)**2)/(1._dp + b*pot0(1:nr)**2)
         pot1(1:nr) = pot1(1:nr)*delta(1:nr)
         dfexc = fourpi*integrate_grid(pot1(1:nr), basis0%grid)
         DEALLOCATE (pot0, pot1)
      CASE DEFAULT
         CPABORT("Unknown functional in atom_dpot_lda")
      END SELECT

      DEALLOCATE (rho0, rho1, drho0, drho1, delta)

      CALL timestop(handle)

   END SUBROUTINE atom_dpot_lda

! **************************************************************************************************
!> \brief Initialise object that holds components needed to compute the exchange-correlation
!>        functional on the atomic radial grid.
!> \param rho_set  electron density object to initialise
!> \param nspins   number of spin components
!> \param needs    components needed to compute the exchange-correlation functional
!> \param rho      electron density on the grid
!> \param drho     first derivative of the electron density on the grid
!> \param tau      kinetic energy density on the grid
!> \param lap      second derivative of the electron density on the grid
!> \param na       number of points on the atomic radial grid
!> \par History
!>    * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
!>    * 08.2008 created as calculate_atom_vxc_lda() [Juerg Hutter]
! **************************************************************************************************
   SUBROUTINE fill_rho_set(rho_set, nspins, needs, rho, drho, tau, lap, na)

      TYPE(xc_rho_set_type), INTENT(INOUT)               :: rho_set
      INTEGER, INTENT(IN)                                :: nspins
      TYPE(xc_rho_cflags_type), INTENT(in)               :: needs
      REAL(dp), DIMENSION(:, :), POINTER                 :: rho, drho, tau, lap
      INTEGER, INTENT(IN)                                :: na

      REAL(KIND=dp), PARAMETER                           :: f13 = (1.0_dp/3.0_dp)

      INTEGER                                            :: ia

      SELECT CASE (nspins)
      CASE (1)
         CPASSERT(.NOT. needs%rho_spin)
         CPASSERT(.NOT. needs%drho_spin)
         CPASSERT(.NOT. needs%norm_drho_spin)
         CPASSERT(.NOT. needs%rho_spin_1_3)
         CPASSERT(.NOT. needs%tau_spin)
         CPASSERT(.NOT. needs%drho)
         ! Give rho to 1/3
         IF (needs%rho_1_3) THEN
            DO ia = 1, na
               rho_set%rho_1_3(ia, 1, 1) = MAX(rho(ia, 1), 0.0_dp)**f13
            END DO
            rho_set%owns%rho_1_3 = .TRUE.
            rho_set%has%rho_1_3 = .TRUE.
         END IF
         ! Give the density
         IF (needs%rho) THEN
            DO ia = 1, na
               rho_set%rho(ia, 1, 1) = rho(ia, 1)
            END DO
            rho_set%owns%rho = .TRUE.
            rho_set%has%rho = .TRUE.
         END IF
         ! Give the norm of the gradient of the density
         IF (needs%norm_drho) THEN
            DO ia = 1, na
               rho_set%norm_drho(ia, 1, 1) = drho(ia, 1)
            END DO
            rho_set%owns%norm_drho = .TRUE.
            rho_set%has%norm_drho = .TRUE.
         END IF
      CASE (2)
         CPASSERT(.NOT. needs%drho)
         CPASSERT(.NOT. needs%drho_spin)
         ! Give the total density
         IF (needs%rho) THEN
            DO ia = 1, na
               rho_set%rho(ia, 1, 1) = rho(ia, 1) + rho(ia, 2)
            END DO
            rho_set%owns%rho = .TRUE.
            rho_set%has%rho = .TRUE.
         END IF
         ! Give the norm of the total gradient of the density
         IF (needs%norm_drho) THEN
            DO ia = 1, na
               rho_set%norm_drho(ia, 1, 1) = drho(ia, 1) + drho(ia, 2)
            END DO
            rho_set%owns%norm_drho = .TRUE.
            rho_set%has%norm_drho = .TRUE.
         END IF
         ! Give rho_spin
         IF (needs%rho_spin) THEN
            DO ia = 1, na
               rho_set%rhoa(ia, 1, 1) = rho(ia, 1)
               rho_set%rhob(ia, 1, 1) = rho(ia, 2)
            END DO
            rho_set%owns%rho_spin = .TRUE.
            rho_set%has%rho_spin = .TRUE.
         END IF
         ! Give rho_spin to 1/3
         IF (needs%rho_spin_1_3) THEN
            DO ia = 1, na
               rho_set%rhoa_1_3(ia, 1, 1) = MAX(rho(ia, 1), 0.0_dp)**f13
               rho_set%rhob_1_3(ia, 1, 1) = MAX(rho(ia, 2), 0.0_dp)**f13
            END DO
            rho_set%owns%rho_1_3 = .TRUE.
            rho_set%has%rho_1_3 = .TRUE.
         END IF
         ! Give the norm of the gradient of rhoa and of rhob separatedly
         IF (needs%norm_drho_spin) THEN
            DO ia = 1, na
               rho_set%norm_drhoa(ia, 1, 1) = drho(ia, 1)
               rho_set%norm_drhob(ia, 1, 1) = drho(ia, 2)
            END DO
            rho_set%owns%norm_drho_spin = .TRUE.
            rho_set%has%norm_drho_spin = .TRUE.
         END IF
         !
      END SELECT

      ! tau part
      IF (needs%tau) THEN
         IF (nspins == 2) THEN
            DO ia = 1, na
               rho_set%tau(ia, 1, 1) = tau(ia, 1) + tau(ia, 2)
            END DO
            rho_set%owns%tau = .TRUE.
            rho_set%has%tau = .TRUE.
         ELSE
            DO ia = 1, na
               rho_set%tau(ia, 1, 1) = tau(ia, 1)
            END DO
            rho_set%owns%tau = .TRUE.
            rho_set%has%tau = .TRUE.
         END IF
      END IF
      IF (needs%tau_spin) THEN
         CPASSERT(nspins == 2)
         DO ia = 1, na
            rho_set%tau_a(ia, 1, 1) = tau(ia, 1)
            rho_set%tau_b(ia, 1, 1) = tau(ia, 2)
         END DO
         rho_set%owns%tau_spin = .TRUE.
         rho_set%has%tau_spin = .TRUE.
      END IF

      ! Laplace
      IF (needs%laplace_rho) THEN
         IF (nspins == 2) THEN
            DO ia = 1, na
               rho_set%laplace_rho(ia, 1, 1) = lap(ia, 1) + lap(ia, 2)
            END DO
            rho_set%owns%laplace_rho = .TRUE.
            rho_set%has%laplace_rho = .TRUE.
         ELSE
            DO ia = 1, na
               rho_set%laplace_rho(ia, 1, 1) = lap(ia, 1)
            END DO
            rho_set%owns%laplace_rho = .TRUE.
            rho_set%has%laplace_rho = .TRUE.
         END IF
      END IF
      IF (needs%laplace_rho_spin) THEN
         CPASSERT(nspins == 2)
         DO ia = 1, na
            rho_set%laplace_rhoa(ia, 1, 1) = lap(ia, 1)
            rho_set%laplace_rhob(ia, 1, 1) = lap(ia, 2)
         END DO
         rho_set%owns%laplace_rho_spin = .TRUE.
         rho_set%has%laplace_rho_spin = .TRUE.
      END IF

   END SUBROUTINE fill_rho_set

! **************************************************************************************************
!> \brief Calculate PADE exchange-correlation (xc) energy functional and potential
!>        in spin-restricted case.
!> \param rho  electron density
!> \param exc  XC energy functional
!> \param vxc  XC potential
!> \par History
!>    * 12.2008 created [Juerg Hutter]
! **************************************************************************************************
   PURE SUBROUTINE lda_pade(rho, exc, vxc)

      REAL(dp), DIMENSION(:), INTENT(IN)                 :: rho
      REAL(dp), DIMENSION(:), INTENT(OUT)                :: exc, vxc

      REAL(KIND=dp), PARAMETER :: a0 = 0.4581652932831429E+0_dp, a1 = 0.2217058676663745E+1_dp, &
         a2 = 0.7405551735357053E+0_dp, a3 = 0.1968227878617998E-1_dp, &
         b1 = 1.0000000000000000E+0_dp, b2 = 0.4504130959426697E+1_dp, &
         b3 = 0.1110667363742916E+1_dp, b4 = 0.2359291751427506E-1_dp, f13 = (1.0_dp/3.0_dp), &
         rsfac = 0.6203504908994000166680065_dp

      INTEGER                                            :: i, n
      REAL(KIND=dp)                                      :: depade, dpv, dq, epade, p, q, rs

      n = SIZE(rho)
      exc(1:n) = 0._dp
      vxc(1:n) = 0._dp

      DO i = 1, n
         IF (rho(i) > 1.e-20_dp) THEN
            rs = rsfac*rho(i)**(-f13)
            p = a0 + (a1 + (a2 + a3*rs)*rs)*rs
            q = (b1 + (b2 + (b3 + b4*rs)*rs)*rs)*rs
            epade = -p/q

            dpv = a1 + (2.0_dp*a2 + 3.0_dp*a3*rs)*rs
            dq = b1 + (2.0_dp*b2 + (3.0_dp*b3 + 4.0_dp*b4*rs)*rs)*rs
            depade = f13*rs*(dpv*q - p*dq)/(q*q)

            exc(i) = epade*rho(i)
            vxc(i) = epade + depade
         END IF
      END DO

   END SUBROUTINE lda_pade

! **************************************************************************************************
!> \brief Calculate PADE exchange-correlation (xc) energy functional and potential
!>        in spin-polarized case.
!> \param rhoa  spin-alpha electron density
!> \param rhob  spin-beta electron density
!> \param exc   XC energy functional
!> \param vxca  spin-alpha XC potential
!> \param vxcb  spin-beta XC potential
!> \par History
!>    * 02.2010 created [Juerg Hutter]
! **************************************************************************************************
   PURE SUBROUTINE lsd_pade(rhoa, rhob, exc, vxca, vxcb)

      REAL(dp), DIMENSION(:), INTENT(IN)                 :: rhoa, rhob
      REAL(dp), DIMENSION(:), INTENT(OUT)                :: exc, vxca, vxcb

      REAL(KIND=dp), PARAMETER :: a0 = 0.4581652932831429E+0_dp, a1 = 0.2217058676663745E+1_dp, &
         a2 = 0.7405551735357053E+0_dp, a3 = 0.1968227878617998E-1_dp, &
         b1 = 1.0000000000000000E+0_dp, b2 = 0.4504130959426697E+1_dp, &
         b3 = 0.1110667363742916E+1_dp, b4 = 0.2359291751427506E-1_dp, &
         da0 = 0.119086804055547E+0_dp, da1 = 0.6157402568883345E+0_dp, &
         da2 = 0.1574201515892867E+0_dp, da3 = 0.3532336663397157E-2_dp, &
         db1 = 0.0000000000000000E+0_dp, db2 = 0.2673612973836267E+0_dp, &
         db3 = 0.2052004607777787E+0_dp, db4 = 0.4200005045691381E-2_dp, f13 = (1.0_dp/3.0_dp), &
         f43 = (4.0_dp/3.0_dp)
      REAL(KIND=dp), PARAMETER :: fxfac = 1.923661050931536319759455_dp, &
         rsfac = 0.6203504908994000166680065_dp

      INTEGER                                            :: i, n
      REAL(KIND=dp)                                      :: dc, dpv, dq, dr, dx, fa0, fa1, fa2, fa3, &
                                                            fb1, fb2, fb3, fb4, fx1, fx2, p, q, &
                                                            rhoab, rs, x, xp, xq

! 1/(2^(4/3) - 2)

      n = SIZE(rhoa)
      exc(1:n) = 0._dp
      vxca(1:n) = 0._dp
      vxcb(1:n) = 0._dp

      DO i = 1, n
         rhoab = rhoa(i) + rhob(i)
         IF (rhoab > 1.e-20_dp) THEN
            rs = rsfac*rhoab**(-f13)

            x = (rhoa(i) - rhob(i))/rhoab
            IF (x < -1.0_dp) THEN
               fx1 = 1.0_dp
               fx2 = -f43*fxfac*2.0_dp**f13
            ELSE IF (x > 1.0_dp) THEN
               fx1 = 1.0_dp
               fx2 = f43*fxfac*2.0_dp**f13
            ELSE
               fx1 = ((1.0_dp + x)**f43 + (1.0_dp - x)**f43 - 2.0_dp)*fxfac
               fx2 = ((1.0_dp + x)**f13 - (1.0_dp - x)**f13)*fxfac*f43
            END IF

            fa0 = a0 + fx1*da0
            fa1 = a1 + fx1*da1
            fa2 = a2 + fx1*da2
            fa3 = a3 + fx1*da3
            fb1 = b1 + fx1*db1
            fb2 = b2 + fx1*db2
            fb3 = b3 + fx1*db3
            fb4 = b4 + fx1*db4

            p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
            q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
            dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
            dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
                                      4.0_dp*fb4*rs)*rs)*rs
            xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
            xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs

            dr = (dpv*q - p*dq)/(q*q)
            dx = 2.0_dp*(xp*q - p*xq)/(q*q)*fx2/rhoab
            dc = f13*rs*dr - p/q

            exc(i) = -p/q*rhoab
            vxca(i) = dc - dx*rhob(i)
            vxcb(i) = dc + dx*rhoa(i)
         END IF
      END DO

   END SUBROUTINE lsd_pade

END MODULE atom_xc
